Ground state of two electrons on a sphere 



Pierre-Frangois Loos and Peter M. W. GilQ 
Research School of Chemistry, Australian National University, Canberra ACT 0200, Australia 

(Dated: February 17, 2010) 

We have performed a comprehensive study of the singlet ground state of two electrons on the 
surface of a sphere of radius R. We have used electronic structure models ranging from restricted 
and unrestricted Hartree-Fock theory to explicitly correlated treatments, the last of which lead to 
near-exact wavefunctions and energies for any value of R. M0ller-Plesset energy corrections (up to 
fifth-order) are also considered, as well as the asymptotic solution in the large-/? regime. 

PACS numbers: 31.15.ac, 31.15.ve, 31.15.xp, 31.15.xr, 31.15.xt 

Keywords: electron correlation, cusp condition, Hartree-Fock solution, symmetry-broken solution, perturba- 
tion theory, M0ller-Plesset, asymptotic expansion, configuration interaction, explicitly correlated methods, 
Hylleraas expansion. 



I. INTRODUCTION 

Exactly (or very accurately) solvable models have 
ongoing value and are valuable both for illuminating 
more complicated systems and for testing theoretical ap- 
proaches, such as density functional methods [THS]- One 
such model is the Hooke's law atom (or Harmonium) 
which is composed of two electrons bound to a nucleus by 
a harmonic potential but repelling Coulombically. This 
system was first considered more than 40 years ago by 
Kestner and Sinanoglu [4], but solved analytically in 1989 
by Kais et al. [5] for a particular value of the harmonic 
force constant and, later, for a countably infinite set force 
constants [6]. 

A related system, studied by Alavi and co-workers 
[7H5], consists of two electrons, interacting through a 
Coulomb potential, but confined within a ball of radius 
R. This possesses a number of interesting features, in- 
cluding the formation of a "Wigner molecule" for large 
R 10 . The spontaneous formation of such molecules 
can also occur in quantum dots and is analogous to the 
Wigner crystallization [11] of the uniform electron gas. 

If the two electrons are constrained to remain on the 
surface of the sphere, one obtains a model that Berry and 
co-workers have used [12 - 15 to understand both weakly 
and strongly correlated systems, such as the ground and 
excited states of the helium atom, and also to suggest the 
"alternating" version of Hund's rule [16]. Seidl studied 
this system in the context of density functional theory 
[T7] in order to test the ISI (interaction-strength interpo- 
lation) model [IS]. For this purpose, he derived accurate 
solutions in both the weak interaction limit (the small 
R regime) and the strong interaction limit (the large R 
regime). He also obtained accurate results by numerical 
integration of the Schrodinger equation. 

In this paper, we are interested in the 1 S ground state 
of two electrons on the surface of a sphere of radius R. 
This allows us to restrict our study to the symmetric 
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spatial part of the wavefunction and ignore the spin co- 
ordinates. We have extended Seidl's analysis and per- 
formed an exhaustive study using a range of models. We 
restrict our analysis to the repulsive potential case; the 
strong-attraction limit (attractive potential) is carefully 
examined in Ref. |17j . 

Restricted and unrestricted Hartree-Fock (HF) solu- 
tions are discussed in Section [Hi] and the strengths and 
weaknesses of M0ller-Plesset (MP) perturbation theory 
[19] in Section IV We consider asymptotic solutions 
for large R in Section [VJ and, in Section |VI[ we ex- 



plore several variational schemes including explicitly cor- 
related techniques [2TFfl2~Tl ) that enforce the cusp condi- 
tion (23 [26] . Atomic units are used throughout. 



II. HAMILTONIAN 

The absolute position of the i-th electron is defined 
by its spherical polar angles Oj = (8i,<j>i). The relative 
position of the electrons is conveniently measured by the 
interelectronic angle 8 which they subtend at the origin. 
These coordinates are related by 

cos 6 = cos B\ cos 62 + sin 8\ sin 82 cos (<fii — ^2) (1) 

and we have < u = \r\ — 7*2] < 2R. 
The Hamiltonian is 



H = f + u- 1 



where 



T = T 1 +T 2 = - 



(2) 



(3) 



is the kinetic operator for both electrons and u 1 is the 
Coulomb operator. In terms of 8, the Hamiltonian is 



1 



* = -w(i +COte T8) + 2R CSC 2 Ui 



in which form it becomes clear that the kinetic and po- 
tential parts of H scale with R~ 2 and Br 1 , respectively. 
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III. HARTREE-FOCK APPROXIMATIONS 



Unrestricted Hartree-Fock 



Restricted Hartree-Fock 



In the HF approximation, each electron feels the mean 
field generated by the other electron [57]. The restricted 
Hartree-Fock (RHF) solution 



$ RHF (n 1; o 2 ) = * K ^(0 1 )^ H±lt (0 2 ) 



RHF i 



,RHF / 



(5) 



When R exceeds a critical value, a second, unrestricted 
HF (UHF) solution develops [3HH3"2] in which the two 
electrons tend to localize on opposite sides of the sphere. 
This is analogous to the UHF description of a dissociating 
H 2 molecule [2"T| . 

To obtain this symmetry- broken solution 



places both electrons in an orbital \|/ RHF that is an eigen- 
function of the Fock operator 



$ UH "(0i^ 2 ) = * 
we expand the orbital as 



UHF/ 



UHF 



(tt- 



RHF . 



«2)| 



R 2 dQ.2 



with dQ.2 = sin 62 d$2 d(f>2 ■ 

By definition, the one-electron basis function 



(6) 



* UHF (0i) = X>^(0i) 



(14) 



(15) 



e=o 



*im(«i) 



i? 



(7) 



where the \I>^(0j) = Y e (9 l )/R = Y( (fl t )/R are zonal 
spherical harmonics. The Fock matrix elements in this 
basis are 



Fi 



where Yf_ m is the spherical harmonic of degree £ and order 
m is an eigenfunction of with eigenvalue £(£+1) / (2R 2 ). 
Using the partial-wave expansion [28 



= h(h + 1) 
2i? 2 



«<,A+ E ^3<^ G ^ ( 16 ) 



1 = R^^TP^cost 



(8) 



and the addition theorem 



e=o 



Pdcos 9) = Y ?m(ni) Y^tfh) (9) (h ^ h) = J Y (l (9) F, 2 (9) Y, 3 (9) sin 9 d9 



where the two-electron integrals are 

GfX = It*" 1 ! <M*i) <Mfe)) (17) 

Using the partial-wave expansion (TSl) and the relation 



it is straightforward to show that 
|*oo(tt 2 )| 2 



R 2 dn 2 



R 



(10) 



The orbital ^00(^1) is thus an eigenfunction of Fi with 
the eigenvalue 1/R. Moreover, it follows from the orthog- 
onality of the spherical harmonics that 



(24 + l)(2£ 2 + l)(2^ + l) fh £ 2 £3 
4tt ^0 

(18) 

between the integrals of three spherical harmonics and 
the Wigner 3j-symbols [53], we find 



G 



R 



E 

1=0 



47T 

2£ + l 



(£i£2£)(£ 3 h£) (19) 



hn.O 



(11) 



which ensures the stationarity of the RHF energy with 
respect to the orbitals ^>e m . 

The ground-state RHF energy is thus 



where selection rules [33j restrict the terms in the sum. 
The UHF energy is then given by 



E 



UHF 



=E3 



E 



RHF 



1 

R 



,£(£+l) 
R 2 



■ E CtiCt-,Fiit 3 ( 2 0) 



and the normalized RHF wavefunction is 



RHF 



AttR 2 



(12) 



(13) 



The first term is the kinetic energy and is positive. How- 
ever, for sufficiently large R, this is outweighed by nega- 
tive contributions in the second term and it is these that 
drive the symmetry-breaking process. 



which yields a uniform electron density over the surface 
of the sphere. 



For computational reasons, we truncate the sum in ( 15 ) 
at £ — L but, for all of the radii R considered in this 
study, we found that L — 15 suffices to obtain £? UHF 
with an accuracy of 1CU 12 . 
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TABLE I: RHF, UHF and exact energies for various R. 
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£.UHF 
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FIG. 1: RxE R ™ (dashed), 7?x£ UHF (dotted) and RxE c 
(solid) as a function of R. 

1.1 




As Table |T] and Figure Q] show, the UHF solut ion be- 
comes lower than the RHF one for R > R CTlt ps 1.5 and 
the UHF, not RHF, energy behaves correctly for large R. 
Specifically, it can be shown that 



lim RE BSF = 1 
lim RE vliF = 1/2 

K— >oo 



(21) 
(22) 



The UHF result reflects the Coulomb interaction be- 
tween two electrons localized on opposite sides of the 
sphere |17j . a phenomenon known as Wigner crystalliza- 
tion [lOj [H] . The difference between the UHF and exact 
energies (i.e. the correlation energy) appears to decay as 

o(ir 3 / 2 ). 



IV. EXPANSION FOR SMALL R 

In M0ller-Plesset perturbation theory, the Hamiltonian 
of the system is partitioned as 



H = H Q + V 



(23) 



where Hq the zeroth-order Hamiltonian and V is a per- 
turbation operator and, in our case, we have 

H =f (24) 

V = u- 1 (25) 

The ground-state wavefunction and energy are expanded 

$ = $(°) + $W + $( 2 ) + $( 3 ) + . . . (26) 

E = + + E& + E& + ... (27) 

We will refer to E^ as the nth-order energy and define 
the MPn correlation energy as 



E 



MP) 







rn=1 



Dimensional analysis reveals that 



$ = 
E = 



9>o 

R 2 ' 
£o_ 

R 2 



R 

£l 

R 



,R 



xR A 



5 R* 

3 



E3R + £<±R + £5-??' 



(28) 

(29) 
(30) 



where the <j) n are functions of and the e n are numbers. 
From (12) and (13), we see </>o = 1/4tt, £0 = and £1 = 1. 
The excited eigenfunctions of Hq are given by 

^(n ll n 2 ) = * <imi (n 1 )^ 2m ,(n 2 ) (3i) 

and we can expand the exact wavefunction $ in this basis. 
However, for the S ground state, angular momentum 
theory [33J|M] limits the combinations of £\, £2, mi and 
ui2 that contribute and it is more efficient to expand $ 
in the basis of two-electron functions 



AitR 2 



Pe(cost 



which are eigenfunctions of f with eigenvalues 

£(£ + 1) 



E e 



R? 



(32) 



(33) 



A. First-order wavefunction 

In the intermediate normalization, the first-order wave- 
function is 



R 



00 ($RHF 

E ' 

1=1 



V 



E — Ee 



1 00 1 



AttR j-^ l{l + 1) 



4 



Using the Legendre generating function 

1 



e=o 



(35) 



the sum in ( 34 1 can be found in closed-form, yielding 
1 



$W(0) = 
or, equivalently, 
*W(u) 



47ri? 



2 In 1 + sin ■ 



1 



1 



4irR 



21nfl + 



(36) 



(37) 



and these yield the normalized first-order wavefunction 



MP1 



v /l + (161n2-ll)i? 2 



(38) 



The true ground-state wavefunction must be nodeless. 
However, it is easy to show that the MP1 wavefunction 
possesses a node if R > 1, leading us to anticipate that 
$ MP1 will be a poor wavefunction for large spheres. 



B. Second- and third-order energies 

According to the Wigner 2n + 1 rule [35] , the lst-order 
wavefunction generates the 2nd- and 3rd-order energies. 
The 2nd-order energy, which has previously been found 
by Seidl [T7], is given by 



E {2) EE S 2 = ( $ 



-,RHF 



V 



$(1) 



4£-3 

-0.227 411 278. 



(39) 



where C = In 2. As Table [TT| shows, the MP2 correlation 
energy is an excellent approximation for small R but, 
because it is independent of R, it is poor for large R. 

It is surprising to find that E^ is so much larger than 
the limiting correlation energies [36) of the helium-like 
ions (-0.0467) or Hooke's Law atoms (-0.0497). 

The 3rd-order energy is given by 



E&> = e 3 R = /*« 
and this yields 

s 3 = 8(£ 2 - 5£ + 3) 
= +0.117 736 889. 



(40) 



(41) 



which agrees with Seidl's rough estimate [17]. Table [TT| 
shows that MP3 gives an improvement over MP2 but 
that it, too, eventually breaks down as R increases. 



C. Second-order wavefunction 

To find the 4th- and 5th-order energies, we need the 
2nd-ordcr wavefunction. This is given by 



* (2) w = E 

which yields 

* (2) (*)^ E E 



V -£« 



Eq — Eg 



(42) 



oo £i+e 2 



e 1 ,£ 2 =ie=\e 1 -e 2 \ 
24 + 1 



A 4(4 + 1)4(4 + 1) \o o y 

Using the identity 



Pi 2 (cos ( 



eJ^^<^<</> 



hi 



(l-x)(l + y) 



for x <y, we eventually obtain 

$ (2) (u) = (2£ 2 - 2£ + 5)0o + (2£ - 5)^i(u) 
1 

+ - 

where Li 2 is the dilogarithm function |37j 



(43) 



1 (44) 



7T 

4 











--) 


-2Li 2 f-^-V 


ARJ 


V 2RJ 



(45) 



D. Fourth- and fifth-order energies 

The Wigner 2n + l rule and the closed- form expression 
of $( 2 ) yield the 4th- and 5th-order coefficients 



e 4 = ^(4£ 3 -42£ 2 - 

= -0.050 275 600 . 

£5 = |(5£ 4 -90£ 3 - 

+ (216-80£)C(3) 
= +0.013 957 832. 



96£-45) - 12C(3) 



450£ 2 - 660£ + 252) 



(46) 



(47) 



where £ is the Riemann zeta function. 

The MPn correlation energies for various values of R 
are reported in Table [IT] and illustrated in Figure [2] The 
results show that MP4 and MP5 are very accurate for 
small R and, indeed, the latter is reasonable up to R ~ 1. 

The MP expansion converges for radii R within the 
radius of convergence 



R cvg 



lim 

n— > 00 



£n+l 



(48) 



From our results, it seems that i? cvg > 2, but it is not 
possible to be more precise than this [T71 [T5] . 
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FIG. 2: MPn correlation energies as a function of R. The 
exact correlation energy is shown as the solid curve. 
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V. EXPANSION FOR LARGE R 

A. Harmonic approximation 

For large R (LR), the potential dominates the kinetic 
energy and the electrons tend to localize on opposite sides 
of the sphere. The classical mechanical energy would be 



E e-e = 



1 

2R 



(49) 



but, quantum mechanically, the kinetic energies of the 
electrons cannot vanish and each electron therefore main- 
tains a zero-point oscillation around its equilibrium po- 
sition with an angular frequency lj. Such phenomena 
are ubiquitous in strongly correlated systems, as demon- 
strated by Seidl and his co-workers [T71 H51 l3"gHlU] . 

In this limit, the supplementary angle £ = it — 6 is the 
natural coordinate and the Hamiltonian becomes 



1 / d 2 d 

H ^-w{de +cot ^ 



— sec - 
2R 2 



For small oscillations (£ ~ 0), the Taylor series 
cot£ = r 1 -£/3-£ 3 /45 + ... 
sec(£/2) = 1 + if/8 + 5£ 4 /384 + 61<f /46080 + . 
yield the harmonic-oscillator Hamiltonian 

whose ground-state wavefunction and energy are 
1 



*o(0 = 



E 



2V2irR 7 / 4 

LRO _ c(0) _ J_ 

t 2i? 



exp(-Vi?r/8) 
1 



2i? 3 / 2 



(50) 



(51) 
(52) 



(53) 

(54) 
(55) 



The second term is the zero-point energy associated with 
harmonic oscillations of angular frequency w = l/i? 3//2 
and it appears that this is the leading error in the UHF 
description at large R. 



By analogy with the small- R expansion (30), we would 



like to construct a large- R asymptotic expansion 



= 5l 
R 

where we know 771 



Hi 
R 2 



(56) 



R 3/2 R 2 R 5/2 ■■■ 

r)2 = 1/2. The nth excited state of 
the Hamiltonian ( 53 1 has the wavefunction and energy 

L n (V^/4)$£(£) (57) 



2,^ 



(58) 



where L„ is the Laguerre polynomial of degree n |29j . 
The anharmonic corrections, £^ and £^ 2 \ can be found 
|41) using the perturbation operators 



R??>~di 4 

R 2 45 d£ 
The first-order correction is 



WW = - 



1 5e 4 

2R 384 
1 61 
2^46 080 



(59) 
(60) 



4>r 



p 00 

Jo 

and this yields 773 = —1/8 and therefore 

^LRl 



1 1 

+ 



2i? 2R?/ 2 
The second-order correction is 

£ (2) = E 



1 

8^2 



n=l 



EX - E± 



W^ 2 



(61) 



(62) 



(63) 



but because of the orthogonality and recurrence relations 
of Laguerre polynomials [H] , only the first two terms in 
the sum in (63) are non-zero and one finds 774 = —1/128 
and therefore 



E 



LR2 



1 



1 



1 



1 



2R 2R?I 2 8R 2 128R 5 / 2 



(64) 



From the results in Table [H] and Figure [3j it seems that 
the asymptotic expansion converges toward the exact en- 
ergy and is reasonably accurate for R > 3. 

Through judicious use of the 5th-order truncation of 
p0| ) and the 2nd-order truncation of ( p36| , one can pre- 
dict satisfactory energies over a wide range of R values. 
However, there remains a region (1 < R < 3) where both 
the small- R and large- R solutions are inadequate. 
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FIG. 3: Correlation energies (relative to RHF) from E LR0 
(dashed), E LR1 (dotted) and £ LR2 (dot-dashed), £; MP5 (small 
dash) and £ exact (solid) as a function of R. 

0.00 r 




TABLE III: Convergence of correlation energies with respect 
to the number L of terms in the CI, R12-CI and Hylleraas 
wavef unctions. All results pertain to the sphere with R = 1. 



VI. VARIATIONAL WAVEFUNCTIONS 
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5 


-0 


145 
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-0 


147 
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-0 


147 


218 


935 


941 


10 


-0 


146 


847 


-0 


147 


217 


796 


-0 


147 


218 


935 


944 


15 


-0 


147 


047 


-0 


147 


218 


405 


-0 


147 


218 


935 


944 


20 


-0 


147 


120 


-0 


147 


218 


631 


-0 


147 


218 


935 


944 


25 


-0 


147 


155 


-0 


147 


218 


738 


-0 


147 


218 


935 


944 


30 


-0 


147 


174 


-0 


147 


218 


797 


-0 


147 


218 


935 


944 


35 


-0 


147 


186 


-0 


147 


218 


833 


-0 


147 


218 


935 


944 


40 


-0 


147 


194 


-0 


147 


218 


857 


-0 


147 


218 


935 


944 



A. Configuration interaction 

We begin with a configuration interaction (CI) treat- 
ment wherein the wavefunction is expanded as 



L 

E- 

=0 



(65) 



in the Legendre polynomial basis set ( 32 1 . The resulting 
energy E^ 1 is the lowest eigenvalue of the CI matrix 



<I», 



H 



> = 



R 2 

t=\ei-t 



1) 



47T 

2£+ 1 



(h he) 



(66) 



where (4 i 2 t) is defined by Q. 

The CI energy as the maximum angular momentum L 
increases is reported in Table [TiT} It converges very slowly 
and even L = 40 yields an accuracy of only 10 -4 . The 
reason for this slow convergence 
satisfy the Kato cusp condition - 



- the failure of (651 to 
is well known. 



B. Hylleraas 

The simplest possible wavefunction with a cusp is 

$ Hy = l + 7« (67) 

which has an explicit linear dependence on the interelec- 
tronic distance u. Kato proved |25) that 7 = 1/2 in nor- 
mal singlet states but, because our electrons are confined 
to a sphere, this 7 does not apply (see below). 
Using the partial-wave expansion 



u = R 



E 

£=0 



1 



1 



2£ + 3 21-1 



P e (cos6) (68) 



one finds that the energy is 



£ Hy (7) = 



6 + 37(7 + 4)i? + 8 7 2 i? 2 



2i?(3 + 8 7 i? + 67 2 i? 2 ) 
and minimizing this with respect to 7 yields 

12 



opt 



7 



9-12R+ V81 + 72i? + 48i? 2 
9 + 12R - V81 + 72i? + 48i? 2 
8i? 2 



(69) 



(70) 



(71) 



Correlation energies obtained from ( 71 1 for several values 



of R are reported in Table [TTJ Despite the simplicity of 
the wavefunction, its energies are surprisingly good with 
a maximum deviation of 0.003 for large R and 0.005 for 
small R. As R tends to zero, the correlation energy ap- 
proaches —0.222 222, which is close to the exact value 
—0.227411. However, as R becomes large, one can show 
that £ Hy - l/(1.58i?) which lies between the RHF and 
UHF energies. The Hylleraas wavefunction is thus a use- 
ful alternative to the small- and large- R solutions in the 
problematic intermediate region (1 < R < 3) with errors 
of 0.000, 0.0011 and 0.0021 for R = 1, 2 and 3, respec- 
tively. 



R12-CI 



Using the Hylleraas wavefunction (67) as the reference 



for a CI expansion yields the R12-CI wavefunction 

L 

$£ 12 - CI ((9) = P$ Hy ((9) +^T £ $ £ (6>) (72) 



where 



(73) 



£=1 



is a projection operator that ensures orthogonality be- 
tween the reference wavefunction and the excited deter- 
minants and / is the identity operator. The coupling 
coefficients between two basis functions are the same as 



those for the conventional CI calculation ( 66 ) but with a 
correction for the matrix element 



^$Hy 



H 



+ 7 



$RHF 

£(£ + 



(74) 



involving the ground state and the excited determinants. 
It is no longer possible to optimize 7 in closed form so 
we used the value given by ( 70 ) . 



As Table [TTT] shows, the R12-CI energies converge much 
more rapidly with L than the CI energies and, for exam- 
ple, Ef 12 ~ cl is more accurate than E^q. This illustrates 
the importance of including a term that is linear in u. 
However, although this term enhances the initial conver- 
gence rate, the asymptotic behavior of the CI and R12-CI 
schemes are identical. Therefore, we now investigate the 
effect of including higher-order u terms. 



where m = min(i, j) and a 



min(i,j) 
max(zj') 



Table 



III 



reveals the remarkable convergence of 2?£ oly . 



Using L = 40 and R = 1, for example, we find 



E 



poly _ 
40 — 



0.852 781 065 056 462 665 400 437 966 038 710 



264 283 589 518 406 360 162 484 313 983 (80) 

The convergence is slower for larger values of R, but still 
impressive. For example, using L — 40 and R — 1000, 
the energy is still correct to 49 digits. The ease with 
which we can obtain these Schrodinger eigenvalues can 
be traced to the fact that the polynomial basis efficiently 
models all of the singularities (the lst-order cusp, the 
third-order cusp, etc.) in the exact wavefunction. 

In recent work |17j . Seidl reported correlation ener- 
gies based on his numerical integration of the Schrodinger 
equation from 9 = to 7r using Q and we have included 
these in Table |H| It appears that sonic of his energies for 
small R are slightly inaccurate. 



VII. CONCLUSION 



D. Polynomial 



In terms of the distance u, the Hamiltonian is 
,2 



H 



(. 



u 



1 



du 2 



\AR 2 

and a Kato-like analysis 



/ 3w 
V4R2 



1\ d 
u j du 



1 

it 



(75) 

reveals the cusp condition 



$'(0) 
*(0) 



= 1 



(76) 



which deviates from the normal value of 1/2 [T7] . 

The natural generalization of the Hylleraas wavefunc- 
tion (67) is a polynomial and it is convenient to select 



the orthonormal basis of Jacobi polynomials 

,(i,o) 



Et(u) 



P 



4ttR 2 

and write the wavefunction as 

L 



RJ 



poly 



(77) 



(78) 



The energy £'£ oly is the lowest eigenvalue of the matrix 



H 



(m 2 — l)(am + 5i,j) am 



4i? 2 



R 



(79) 



In this article, we have reported results for the ground 
state of a simple two-electron system that is described 
by a single parameter R. Although we cannot solve its 
Schrodinger equation in closed form, we have found accu- 
rate wavefunctions and energies for small R (the weakly 
correlated limit) and large R, (the strongly correlated 
limit). For R -C 1, M0ller-Plesset perturbation theory 
yields results close to the exact solution; for R 3> 1, ac- 
curate results can be found by considering the zero-point 
oscillations of the appropriate Wigner molecule. 

We have also explored variational schemes that yield 
satisfactory results for all R. In particular, we have dis- 
covered a polynomial wavefunction that easily yields re- 
sults of any required accuracy. 

We believe that our results will be useful in the future 
development of accurate correlation functionals within 
density-functional theory [321 ESj aud intracule functional 
theory 
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